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Abstract 

We present the first calculations of the gravitational radiation produced by 

nonaxisymmetric dynamical instability in a rapidly rotating compact star. 

The star deforms into a bar shape, shedding ~ 4% of its mass and ~ 17% 

of its angular momentum. The gravitational radiation is calculated in the 

quadrupole approximation. For a mass M ~ 1.4 Mq and radius R ~ 10 km, 

the gravitational waves have frequency ~ 4 kHz and amplitude h ~ 2 x 10 -22 

at the distance of the Virgo Cluster. They carry off energy AE/M ~ 0.1% 

and radiate angular momentum AJ/J ~ 0.7%. 
PACS numbers: 04.30,+x, 04.80. +z, 97.60.-s 
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The detection of gravitational waves from astrophysical sources is one of the key frontiers 
for research in general relativity [|I|] . Among the sources that might be observed by detectors 
currently under development are those due to global rotational instabilities in collapsing or 
compact stars For example, a rapidly rotating stellar core that has exhausted its 

nuclear fuel and is prevented from collapsing to neutron star size by centrifugal forces could 
become unstable, potentially shedding enough angular momentum to allow it to collapse 
and form a supernova [§]]. In addition, neutron stars spun up by accretion of mass from 
a binary companion could possibly reach fast enough rotation rates to go unstable |jU§. 
The prospect of several gravitational wave detectors becoming operational within a decade 
means that detailed modeling of such sources and the radiation they produce has a high 
priority. 

Global rotational instabilities arise in fluids from growing modes e ±tm< ^, where m = 2 
is the so-called "bar mode". They are conveniently parametrized by (3 = T/IW 7 !, where T 
is the rotational kinetic energy and W is the gravitational potential energy; see for 
reviews. We focus on the bar instability since it is expected to be the fastest growing mode. 
There are two different physical mechanisms by which this instability can develop. The 
dynamical bar instability is driven by Newtonian hydrodynamics and gravity. It operates 
for fairly large values (3 > (3d and develops on the timescale of about a rotation period. The 
secular instability is due to dissipative processes such as gravitational radiation reaction and 
occurs for typically lower values (3 S < (3 < (3d- It develops on a timescale of several rotation 
periods or longer || . For the uniformly rotating, incompressible, constant density Maclaurin 
spheroids (3 S ~ 0.14 and [3d ~ 0.27. Differentially rotating fluids with a polytropic equation 
of state P oc p 1+1 / n , where n is the polytropic index, are believed to undergo secular and 
dynamical bar instabilities at about these same values of (3 |7]|| . 

Astrophysical sources driven by rotational instabilities are nonlinear, time-dependent, 
and fully 3-dimensional systems; calculation of the gravitational radiation they produce re- 
quires numerical modeling. We have carried out computer simulations of a differentially 
rotating compact star with a polytropic equation of state undergoing the dynamical bar 
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instability. This instability has previously been modeled numerically by Tohline and col- 
laborators in the context of star formation |7|flOHT2| . Our work is the first to calculate the 
gravitational radiation produced by this instability, including waveforms and luminosities. 
It is also a significant advance over the earlier studies because, in addition to using bet- 
ter numerical techniques, we model the fluid correctly using an energy equation. This is 
essential due to the generation of entropy by shocks during the later stages of the evolution. 

The gravitational radiation is calculated using Newtonian gravity in the quadrupole 
approximation [13|. The gravitational waveforms are given by the transverse-traceless (TT) 



components of the metric perturbation r hj? = 2 J^- TT , where = / p(xiXj — \5ijT 2 ) d 3 r 
is the trace-reduced quadrupole moment of the source and we use units in which c—G — 1. 
Here, r 2 = x 2 + y 2 + z 2 , spatial indices i,j = 1,2,3, and a dot indicates a time derivative 
dj dt. For an observer located on the axis at 9 = 0, = in spherical coordinates centered on 
the source, the waveforms take the simple form rh TT = (I xx — f yy )e+ + 2I xy e x , where e + and 
e x are basis tensors for the two polarization states The gravitational wave luminosity 
is L = / an d the angular momentum lost is dJi/dt = -e^ /ij^fj^V where the 
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superscript (3) indicates 3 time derivatives, there is an implied sum on repeated indices, and 
the angle-brackets indicate an average over several wave periods. For these burst sources 
such averaging is not well-defined; therefore we display the unaveraged quantities |f®f® 
and |ejjA:-fjm^im below. The energy emitted as gravitational radiation is AE = J L dt and 
the angular momentum carried away by the waves is AJj = J(dJi/dt)dt. 

The equations of hydrodynamics are integrated using the techniques of smoothed particle 
hydrodynamics (SPH) ||15|| , in which the fluid is modeled as a collection of particles having 
non-zero extent given by a smoothing length h. The value of any physical field is then 
obtained by averaging over all particles within 2h of a given point using kernel estimation. 
We have used the implementation of SPH by Hernquist and Katz known as TREESPH 
I6| , which has variable smoothing lengths and individual particle timesteps. This code has 



been well-tested by its developers; see |TB| for details. To reduce the numerical noise that 
can arise from taking derivatives of iy numerically, we calculate fj,- analytically from the 



SPH equations of motion and obtain the gravitational waveforms directly from quantities 
already available in the code (cf. |I7||); this produces very smooth profiles. The gravitational 
radiation calculation and the use of artificial viscosity to handle shocks have been extensively 
tested and results are available in []16| and ||1~8 |. 



Our calculations begin with an initially axisymmetric differentially rotating star with a 
polytropic equation of state. This equilibrium model is produced using the self-consistent 
field method of Smith and Centrella fl9|], which is based on earlier work of Ostriker and 
Mark f2(| and Hachisu . This method is an iterative technique carried out in cylindrical 



coordinates (w,z,4>). A rotation law is specified by giving the specific angular momentum 
as a function of mass interior to a cylinder of radius w. Following the convention of earlier 
work [|10| - |i2| , |22|| , we use the rotation law for the rigidly rotating, uniform density Maclaurin 



spheroids. Since polytropes do not have uniform density, this produces differentially rotating 
models. The resulting density p(w, z) and angular velocity Q(w) are then converted into a 
particle model to be evolved with TREESPH. 

The star has polytropic index n = 3/2, total mass M, equatorial radius R, and j3 = 
0.30. This rapidly rotating model is highly flattened, with polar radius ~ 0.21R. Time is 
measured in units of the dynamical timescale for a sphere of radius R, tjj = (R 3 / M) 1//2 . No 
perturbations are imposed upon the initial axisymmetric rotating equilibrium model; the 
dynamical bar instability grows spontaneously from small deviations from axisymmetry in 
the particle model. We carried out runs with N = 2000,4000,8000, and 16,000 particles. 
All models were run for 20to, conserved total energy to ~ 1% and angular momentum 
to ~ 0.1%, and produced generally similar results. Table | shows how some of the bulk 
properties of the star vary with particle number. We ran a set of tests to investigate the 
effect of artificial viscosity on the growth of the bar instability, comparing our results with the 
linearized tensor virial analysis fllQfl . Since the artificial viscosity used in ref. |18| produced 
results that most closely matched the analytic ones, we used this form in the runs discussed 
here. Both the growth rate and rotation speed of the m = 2 bar mode in all our runs were 
closer to the analytic results than the previous numerical studies |T0 |. 



Fig. [i] shows the particle positions for the case N = 16, 000 projected onto the x — y plane 
at 6 different times during the calculation. The star is rotating counterclockwise about the 
z axis. The simulation starts at t — 0, Fig. |l](a). By t = 6.6to, Fig. |l|(b), nonaxisymmetric 
structure is apparent. As this bar shaped structure develops, the amplitude of the m = 2 
mode grows exponentially until t ~ lOtp. During this time a spiral arm pattern develops, 
with mass being shed from the ends of the bar, as seen in Fig. |l](c) and |l|(d). The bar 
and spiral arms exert gravitational torques that cause angular momentum to be transported 
outward [SUM! . The spiral arms expand supersonically and merge together, causing shock 
heating and dissipation; see Fig. |l](e). The system evolves toward a nearly axisymmetric 
state with a core of equatorial radius ~ R and (3 ~ 0.26, as shown in Fig. |l|(f); c.f. p3 



The halo contains ~ 4% of the total mass and ~ 17% of the angular momentum, and is 
bound. Throughout its evolution the entire system remains flattened, with final polar radius 
~ 0.25R. 

The gravitational waveforms rh + (solid line) and rh x (dashed line) for this run are shown 
in Fig. [2] for an observer on the axis at 9 = 0, = 0. Fig. |3] shows the gravitational wave 
luminosity L and Fig. ^ shows the energy AE/M emitted as gravitational waves. The rate 
dJ z /dt at which angular momentum is carried away by the gravitational waves is given in 
Fig. [| and the angular momentum lost to gravitational radiation A J z j J is shown in Fig. |6|. 
Scaling these results to a mass M ~ 1.4 M and radius R ~ 10 km, we find that the 
gravitational radiation has frequency / ~ 4 kHz. The dimensionless amplitude of the waves 
is h ~ 4 x 10~ 19 at a distance r ~ 10 kpc, typical of sources within the Milky Way, and 
h ~ 2 x 10~ 22 at a distance of r ~ 20 Mpc, typical of sources in the Virgo Cluster. The 
energy lost to gravitational radiation is AE/M ~ 0.1% and the angular momentum radiated 
away is A J/ J ~ 0.7%. 

In this Letter we present the first calculations of the gravitational radiation produced 
by the dynamical bar instability. This is an important step in understanding astrophysi- 
cal sources driven by global rotational instabilities. The dynamical instability may be the 
mechanism by which a rapidly rotating, contracting stellar core sheds enough angular mo- 



mentum to allow collapse to neutron star densities |[24|| , and methods developed in this work 
can be applied to study such collapses. Neutron stars spun up by accretion are expected 
to be subject to the secular instability ||. Modeling this process will require including the 
effects of gravitational radiation reaction [^5j and viscous fluids [^6] in the hydrodynamical 
equations. 

We thank L. Hernquist for supplying a copy of TREESPH and S. McMillan, K. Thorne, 
and C. Cutler for interesting and helpful discussions. This work was supported by NSF 
grant PHY-9208914 and the calculations were performed at the Pittsburgh Supercomputing 
Center. 
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FIGURES 

FIG. 1. Particle positions are shown projected onto the x — y plane for various times in the 
evolution of the model with N = 16, 000 particles. 

FIG. 2. Gravitational waveforms rh + (solid line) and rh x (dashed line) are shown for the run 
with N = 16, 000. 

FIG. 3. The gravitational wave luminosity L (M/i?) -5 is shown. This profile has been 
smoothed using simple averaging over a fixed interval of O.Ud centered on each point. 

FIG. 4. The energy [AE/M] (M / R)~ 7 1 2 emitted in the form of gravitational waves is given. 

FIG. 5. The angular momentum [dJ z /dt] M^ 1 (M/R)' 7 ^ 2 carried away by gravitational waves 
is shown. 

FIG. 6. The total angular momentum [AJ Z /J] (M/i?) -5 / 2 lost to gravitational radiation is 
shown. 
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TABLES 

TABLE I. The mass and angular momentum shed to the halo, and the energy and angular 
momentum carried by the waves (for M = 1.4 Mq and R = 10 km), for runs with different particle 



number. 

N 2000 4000 8000 16,000 

mass in halo 6.0% 6.1% 5.7% 4.4% 

J in halo 19% 19% 18% 17% 

AE/M .057% .075% .075% .10% 

AJ Z /J .36% .52% .52% .72% 
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